library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)library(readr)library(tidyr)# Read the CAPTAIN2 EDGE2 RDS fileCAPTAIN2_EDGE2_data <-readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))# Analyze non-zero cells in the prioritization outputcat("Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:\n")
Analyzing cells with non-zero priority values in CAPTAIN2 EDGE2 output:
cat("Cells with zero priority:", zero_cells, " (", round(zero_cells/total_cells*100, 2), "%)\n", sep="")
Cells with zero priority:226942 (97.58%)
Code
cat("Cells with NA priority:", na_cells, " (", round(na_cells/total_cells*100, 2), "%)\n", sep="")
Cells with NA priority:0 (0%)
Code
# Summary statistics of priority valuespriority_summary <-summary(CAPTAIN2_EDGE2_data$Priority)cat("\nSummary statistics of priority values:\n")
Summary statistics of priority values:
Code
print(priority_summary)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00000 0.00000 0.00000 0.02331 0.00000 1.00000
Code
# Distribution of non-zero priority valuesnonzero_priority <- CAPTAIN2_EDGE2_data$Priority[CAPTAIN2_EDGE2_data$Priority >0]cat("\nDistribution of non-zero priority values:\n")
# Load one of your input raster files to extract the correct grid structureraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 data based on PUIDCAPTAIN2_EDGE2_data_with_coords <- CAPTAIN2_EDGE2_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_EDGE2_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_EDGE2_data_nonzero <- CAPTAIN2_EDGE2_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_EDGE2_sf <-st_as_sf( CAPTAIN2_EDGE2_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_EDGE2 <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_EDGE2 <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_EDGE2 <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_EDGE2_plot <-ggplot() +geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_EDGE2# Display the plotprint(CAPTAIN2_EDGE2_plot)
Code
# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_2ndrun_EDGE2_priorities_01_2.png"),plot = CAPTAIN2_EDGE2_plot,width =10,height =6,dpi =300,bg ="white")
FUSE continental 0.1 budget
Code
library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read the CAPTAIN2 FUSE RDS fileCAPTAIN2_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 FUSE data based on PUIDCAPTAIN2_FUSE_data_with_coords <- CAPTAIN2_FUSE_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_FUSE_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_FUSE_data_nonzero <- CAPTAIN2_FUSE_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_FUSE_sf <-st_as_sf( CAPTAIN2_FUSE_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_FUSE <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_FUSE <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_FUSE <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_FUSE_plot <-ggplot() +geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_FUSE, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_FUSE, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_FUSE# Display the plotprint(CAPTAIN2_FUSE_plot)
Code
# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_FUSE_priorities_01_2.png"),plot = CAPTAIN2_FUSE_plot,width =10,height =6,dpi =300,bg ="white")
IUCN continental 0.1 budget
Code
library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read the CAPTAIN2 IUCN RDS fileCAPTAIN2_IUCN_data <-readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 IUCN data based on PUIDCAPTAIN2_IUCN_data_with_coords <- CAPTAIN2_IUCN_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_IUCN_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 IUCN data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_IUCN_data_nonzero <- CAPTAIN2_IUCN_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_IUCN_sf <-st_as_sf( CAPTAIN2_IUCN_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_IUCN <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_IUCN <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_IUCN <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Create the plotCAPTAIN2_IUCN_plot <-ggplot() +geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_IUCN, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_IUCN, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_IUCN# Display the plotprint(CAPTAIN2_IUCN_plot)
Code
# Save the plotggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_priorities_01_2.png"),plot = CAPTAIN2_IUCN_plot,width =10,height =6,dpi =300,bg ="white")
Difference maps
Code
library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)# Read all three index RDS filesCAPTAIN2_IUCN_data <-readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))CAPTAIN2_EDGE2_data <-readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))CAPTAIN2_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))# Load one of your input raster files to extract the correct grid structureraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Create a grid of coordinates for each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Join all three datasets with coordinatesIUCN_with_coords <- CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority) %>%left_join(coords, by ="PUID")EDGE2_with_coords <- CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority) %>%left_join(coords, by ="PUID") FUSE_with_coords <- CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority) %>%left_join(coords, by ="PUID")# Combine all datasetsall_indices <- coords %>%left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority), by ="PUID") %>%left_join(CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority), by ="PUID") %>%left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority), by ="PUID")# Calculate differencesall_indices <- all_indices %>%mutate(# Replace NA with 0 for calculation purposesIUCN =ifelse(is.na(IUCN), 0, IUCN),EDGE2 =ifelse(is.na(EDGE2), 0, EDGE2),FUSE =ifelse(is.na(FUSE), 0, FUSE),# Calculate differencesIUCN_minus_FUSE = IUCN - FUSE,IUCN_minus_EDGE2 = IUCN - EDGE2,EDGE2_minus_FUSE = EDGE2 - FUSE )# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected <-st_transform(world, crs = mcbryde_thomas_2)# Create base thememy_theme <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Filter to non-zero differences for each comparison to reduce plot size# IUCN - FUSEIUCN_FUSE_diff <- all_indices %>%filter(IUCN_minus_FUSE !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# IUCN - EDGE2IUCN_EDGE2_diff <- all_indices %>%filter(IUCN_minus_EDGE2 !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# EDGE2 - FUSEEDGE2_FUSE_diff <- all_indices %>%filter(EDGE2_minus_FUSE !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# Create a diverging color palette for difference maps# Blue for negative (first index lower), white for zero, red for positive (first index higher)diff_colors <-c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0", "#FFFFFF", "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")# 1. IUCN - FUSE Difference MapIUCN_FUSE_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(IUCN - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="IUCN Index minus FUSE Index",x =NULL, y =NULL) + my_theme# 2. IUCN - EDGE2 Difference MapIUCN_EDGE2_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(IUCN - EDGE2)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="IUCN Index minus EDGE2 Index",x =NULL, y =NULL) + my_theme# 3. EDGE2 - FUSE Difference MapEDGE2_FUSE_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(EDGE2 - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="EDGE2 Index minus FUSE Index",x =NULL, y =NULL) + my_theme# Display all plotsprint(IUCN_FUSE_plot)
Code
print(IUCN_EDGE2_plot)
Code
print(EDGE2_FUSE_plot)
Code
# Save all plotsggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_minus_FUSE_difference_2.png"),plot = IUCN_FUSE_plot,width =10,height =6,dpi =300,bg ="white")ggsave(filename = here::here("outputs", "CAPTAIN2_IUCN_minus_EDGE2_difference_2.png"),plot = IUCN_EDGE2_plot,width =10,height =6,dpi =300,bg ="white")ggsave(filename = here::here("outputs", "CAPTAIN2_EDGE2_minus_FUSE_difference_2.png"),plot = EDGE2_FUSE_plot,width =10,height =6,dpi =300,bg ="white")# Optionally, create a panel with all three difference mapslibrary(patchwork)# Combine all plotsall_diffs_plot <- IUCN_FUSE_plot / IUCN_EDGE2_plot / EDGE2_FUSE_plot +plot_annotation(title ="Differences Between Conservation Priority Indices",subtitle ="Budget: 0.1, Replicates: 100",theme =theme(plot.title =element_text(hjust =0.5),plot.subtitle =element_text(hjust =0.5)) )#all_diffs_plot# Save the combined plot#ggsave(# filename = here::here("outputs", "CAPTAIN2_all_differences.png"),# plot = all_diffs_plot,# width = 10,# height = 15,# dpi = 300,# bg = "white"#)
Species level priorities
Code
# Load required packageslibrary(here)library(readr)library(dplyr)library(tidyr)library(ggplot2)# Read the protected range fractions RDS fileprotected_fractions <-readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions_2ndrun.rds"))# Read the continental shark conservation metrics CSV fileshark_metrics <-read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))# Rename column in shark_metrics to match bettershark_metrics <- shark_metrics %>%rename(Species =`Species name`)# Order shark_metrics alphabetically by Species nameshark_metrics <- shark_metrics %>%arrange(Species)# Add an original order ID to protected_fractions to maintain its original orderprotected_fractions$original_order <-1:nrow(protected_fractions)# Add row number as IDs to both datasetsprotected_fractions$protected_ID <-1:nrow(protected_fractions)shark_metrics$species_ID <-1:nrow(shark_metrics)# Check if the datasets have the same number of rowsif(nrow(protected_fractions) ==nrow(shark_metrics)) {# Create index mapping - this maintains the original protection data ordering# while allowing us to associate with alphabetically ordered species names indices <-data.frame(protected_ID =1:nrow(protected_fractions),species_ID =1:nrow(shark_metrics) )# Join protected_fractions with indices protected_with_indices <- protected_fractions %>%left_join(indices, by ="protected_ID")# Join shark_metrics with indices shark_with_indices <- shark_metrics %>%left_join(indices, by ="species_ID")# Now join the datasets, matching on species_ID and protected_ID combined_data <- protected_with_indices %>%inner_join( shark_with_indices,by =c("species_ID", "protected_ID"),suffix =c("_captain", "_original") ) %>%# Sort by the original order of protected_fractions arrange(original_order)cat("Successfully joined datasets with", nrow(combined_data), "species\n")cat("First few species in combined dataset:\n")print(head(combined_data[, c("Species_captain", "Species_original")]))# Define IUCN categories and order - using only the first 5 categories iucn_labels <-c("1"="LC", "2"="NT", "3"="VU", "4"="EN", "5"="CR" ) iucn_order <-c("LC", "NT", "VU", "EN", "CR")# Define colors for IUCN categories iucn_colors <-c("LC"="#50C878", # Green"NT"="#FFFF00", # Yellow"VU"="#FFA500", # Orange"EN"="#FF8C00", # Dark Orange"CR"="#FF0000"# Red )# 1. IUCN Boxplot iucn_boxplot <- combined_data %>%mutate(IUCN_status =factor(IUCN_original, levels =1:5, labels = iucn_order),protection_percentage = IUCN_captain *100 ) %>%ggplot(aes(x = IUCN_status, y = protection_percentage)) +# geom_violin(aes(fill = IUCN_status, color = IUCN_status), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="IUCN Priority Index: Range Protection by IUCN Status",x ="IUCN Red List threat status", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_fill_manual(values = iucn_colors) +scale_color_manual(values = iucn_colors) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# 2. FUSE Boxplot fuse_boxplot <- combined_data %>%mutate(FUSE_category =factor(FUSE_original),protection_percentage = FUSE_captain *100 ) %>%ggplot(aes(x = FUSE_category, y = protection_percentage)) +# geom_violin(aes(fill = FUSE_category, color = FUSE_category), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="FUSE Priority Index: Range Protection by FUSE Score",x ="FUSE Score", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# 3. EDGE2 Boxplot edge2_boxplot <- combined_data %>%mutate(EDGE2_category =factor(EDGE2_original),protection_percentage = EDGE2_captain *100 ) %>%ggplot(aes(x = EDGE2_category, y = protection_percentage)) +# geom_violin(aes(fill = EDGE2_category, color = EDGE2_category), # trim = FALSE, # alpha = 0.5) +geom_jitter(width =0.1, size =0.6, alpha =0.5, color ="darkgray") +geom_boxplot(width =0.1, fill ="white", color ="black", outlier.shape =NA, alpha =0.8) +labs(title ="EDGE2 Priority Index: Range Protection by EDGE2 Score",x ="EDGE2 Score", y ="Range protected (%)") +theme_minimal() +theme(plot.title =element_text(hjust =0.5),legend.position ="none",panel.grid.major.x =element_blank() ) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# Display individual plotsprint(iucn_boxplot)print(fuse_boxplot)print(edge2_boxplot)# Save individual plotsggsave(here::here("outputs", "iucn_protection_boxplot.png"), iucn_boxplot, width =10, height =6, dpi =300)ggsave(here::here("outputs", "fuse_protection_boxplot.png"), fuse_boxplot, width =10, height =6, dpi =300)ggsave(here::here("outputs", "edge2_protection_boxplot.png"), edge2_boxplot, width =10, height =6, dpi =300)# Create a nicely formatted table showing the ordered species species_table <- combined_data %>% dplyr::select(Species_captain, Species_original, IUCN_captain, FUSE_captain, EDGE2_captain, IUCN_original, FUSE_original, EDGE2_original) %>%arrange(Species_original) %>%head(20) # Just show the first 20 for display# Print species tablecat("\nFirst 20 species (alphabetically by original species name):\n")print(species_table)# Save the full combined datawrite.csv(combined_data, here::here("outputs", "combined_species_data.csv"), row.names =FALSE)} else {cat("ERROR: Datasets have different number of rows.\n")cat("Protected fractions:", nrow(protected_fractions), "rows\n")cat("Shark metrics:", nrow(shark_metrics), "rows\n")}
Successfully joined datasets with 1000 species
First few species in combined dataset:
Species_captain Species_original
1 Species_1 Acroteriobatus_annulatus
2 Species_2 Acroteriobatus_blochii
3 Species_3 Acroteriobatus_leucospilus
4 Species_4 Acroteriobatus_ocellatus
5 Species_5 Acroteriobatus_salalah
6 Species_6 Acroteriobatus_variegatus
# Load required librarieslibrary(patchwork)library(ggplot2)# Create a function to add labels (A), (B), etc.add_panel_labels <-function(plot, label) { plot +theme(plot.title =element_text(face ="bold", hjust =0, size =12) ) +labs(title =paste0("(", label, ")"))}# Add labels to each plot# First gridCAPTAIN2_IUCN_msmap_labeled <-add_panel_labels(CAPTAIN2_IUCN_msmap, "A")CAPTAIN2_FUSE_msmap_labeled <-add_panel_labels(CAPTAIN2_FUSE_msmap, "B")CAPTAIN2_EDGE2_msmap_labeled <-add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")# Combine plots into two separate grids, each with 3 rows and 1 columngrid1 <- CAPTAIN2_IUCN_msmap_labeled / CAPTAIN2_FUSE_msmap_labeled / CAPTAIN2_EDGE2_msmap_labeled# Display each grid separatelygrid1
Code
# Save the grids if neededggsave(here::here("grid1_maps_continental_2ndrun.png"), grid1, width =8, height =15, dpi =300)
Difference maps
Code
# 1. IUCN - FUSE Difference MapIUCN_FUSE_msmap <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(IUCN - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Difference in Conservation Priorities",#subtitle = "IUCN Index minus FUSE Index",x =NULL, y =NULL) + my_theme# 2. IUCN - EDGE2 Difference MapIUCN_EDGE2_msmap <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(IUCN - EDGE2)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Difference in Conservation Priorities",#subtitle = "IUCN Index minus EDGE2 Index",x =NULL, y =NULL) + my_theme# 3. EDGE2 - FUSE Difference MapEDGE2_FUSE_msmap <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority\n(EDGE2 - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Difference in Conservation Priorities",#subtitle = "EDGE2 Index minus FUSE Index",x =NULL, y =NULL) + my_theme# Load required librarieslibrary(patchwork)library(ggplot2)# Create a function to add labels (A), (B), etc.add_panel_labels <-function(plot, label) { plot +theme(plot.title =element_text(face ="bold", hjust =0, size =12) ) +labs(title =paste0("(", label, ")"))}# Add labels to each plot# First gridIUCN_FUSE_msmap_labeled <-add_panel_labels(IUCN_FUSE_msmap, "A")IUCN_EDGE2_msmap_labeled <-add_panel_labels(IUCN_EDGE2_msmap, "B")EDGE2_FUSE_msmap_labeled <-add_panel_labels(EDGE2_FUSE_msmap, "C")# Combine plots into two separate grids, each with 3 rows and 1 columngrid2 <- IUCN_FUSE_msmap_labeled / IUCN_EDGE2_msmap_labeled / EDGE2_FUSE_msmap_labeled# Display each grid separatelygrid2
Code
# Save the grids if neededggsave(here::here("grid2_maps_continental_2ndrun.png"), grid2, width =8, height =15, dpi =300)
Congruence maps
Code
# Method 1: If you have separate dataframes for each index# Assuming your data has columns: PUID, Priority, and geometry# First, identify high priority areas (>0.9) for each indexhigh_priority_EDGE2 <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$Priority >0.9, ]high_priority_FUSE <- CAPTAIN2_FUSE_sf[CAPTAIN2_FUSE_sf$Priority >0.9, ]high_priority_IUCN <- CAPTAIN2_IUCN_sf[CAPTAIN2_IUCN_sf$Priority >0.9, ]# Find congruent areas (present in all three)# Method using PUID (assuming you have PUID column)congruent_PUIDs <-intersect(intersect(high_priority_EDGE2$PUID, high_priority_FUSE$PUID), high_priority_IUCN$PUID)# Extract congruent areascongruent_areas <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$PUID %in% congruent_PUIDs, ]# Print summarycat("High priority areas (>0.9):\n")
cat("Percentage of overlap:", round(nrow(congruent_areas)/min(nrow(high_priority_EDGE2), nrow(high_priority_FUSE), nrow(high_priority_IUCN))*100, 2), "%\n")
Percentage of overlap: 33.19 %
Code
# Method 2: If you need to merge dataframes first# Create a combined dataframe with all three priority scores# First, extract just the data (without geometry) from the other sf objectsFUSE_data <-st_drop_geometry(CAPTAIN2_FUSE_sf[, c("PUID", "Priority")])IUCN_data <-st_drop_geometry(CAPTAIN2_IUCN_sf[, c("PUID", "Priority")])# Merge with the EDGE2 sf object (keeping geometry)combined_priorities <-merge(CAPTAIN2_EDGE2_sf, FUSE_data, by ="PUID", suffixes =c("_EDGE2", "_FUSE"))combined_priorities <-merge(combined_priorities, IUCN_data, by ="PUID")names(combined_priorities)[names(combined_priorities) =="Priority"] <-"Priority_IUCN"# Identify congruent high priority areascongruent_areas_v2 <- combined_priorities[combined_priorities$Priority_EDGE2 >0.9& combined_priorities$Priority_FUSE >0.9& combined_priorities$Priority_IUCN >0.9, ]# Create a map showing only congruent areascongruent_map <-ggplot() +geom_sf(data = congruent_areas, aes(fill ="Congruent High Priority"), color ="red", size =0.8, alpha =0.8) +geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill =NA, color ="black", size =0.5) +scale_fill_manual(values =c("Congruent High Priority"="red"),name ="Priority Areas") +labs(title ="Congruent High Priority Areas",subtitle ="Areas with Priority > 0.9 in EDGE2, FUSE, and IUCN indices",x =NULL, y =NULL) + my_theme_CAPTAIN2_EDGE2# Display the mapprint(congruent_map)
Code
# Create detailed congruence categories# Define logical conditions for each indexEDGE2_high <- combined_priorities$Priority_EDGE2 >0.9FUSE_high <- combined_priorities$Priority_FUSE >0.9IUCN_high <- combined_priorities$Priority_IUCN >0.9# Create detailed congruence categoriescombined_priorities$congruence_category <-case_when( EDGE2_high & FUSE_high & IUCN_high ~"All three indices", EDGE2_high & FUSE_high &!IUCN_high ~"EDGE2 + FUSE", EDGE2_high &!FUSE_high & IUCN_high ~"EDGE2 + IUCN", !EDGE2_high & FUSE_high & IUCN_high ~"FUSE + IUCN", EDGE2_high &!FUSE_high &!IUCN_high ~"EDGE2 only",!EDGE2_high & FUSE_high &!IUCN_high ~"FUSE only",!EDGE2_high &!FUSE_high & IUCN_high ~"IUCN only",TRUE~"No high priority")# Convert to factor with desired ordercombined_priorities$congruence_category <-factor(combined_priorities$congruence_category,levels =c("All three indices", "EDGE2 + FUSE", "EDGE2 + IUCN", "FUSE + IUCN","EDGE2 only", "FUSE only", "IUCN only"))# Filter data to only include high priority areashigh_priority_data <- combined_priorities[combined_priorities$congruence_category %in%c("All three indices", "EDGE2 + FUSE", "EDGE2 + IUCN", "FUSE + IUCN", "EDGE2 only", "FUSE only", "IUCN only"), ]# Debug: Check the datacat("Data check:\n")
Data check:
Code
cat("Total high priority areas:", nrow(high_priority_data), "\n")
All three indices EDGE2 + FUSE EDGE2 + IUCN FUSE + IUCN
1459 275 1130 75
EDGE2 only FUSE only IUCN only <NA>
154 17 29 6
Code
# Calculate percentagestotal_high_priority <-sum(congruence_summary[names(congruence_summary) !="No high priority"], na.rm =TRUE)congruence_percentages <-round(congruence_summary / total_high_priority *100, 2)cat("\nPercentages of high priority areas:\n")
Percentages of high priority areas:
Code
print(congruence_percentages[names(congruence_percentages) !="No high priority"])
All three indices EDGE2 + FUSE EDGE2 + IUCN FUSE + IUCN
46.48 8.76 36.00 2.39
EDGE2 only FUSE only IUCN only <NA>
4.91 0.54 0.92
Relationship with fishing effort : bar plots per country
Code
library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(raster)library(tidyr)library(ggrepel)library(ggsci) # For nice color palettes# Read all priority dataCAPTAIN2_EDGE2_data <-readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))CAPTAIN2_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))CAPTAIN2_IUCN_data <-readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))# Load fishing dataload(here::here("Data", "Raw", "Predicted_Fishing_Hours_05Deg.Rdata"))# Load marine ecoregion shapefilemeow_ecos <-st_read(here("Data", "Shapefiles", "meow_ecos", "meow_ecos.shp"), quiet =TRUE)# Load raster to get coordinate systemraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")r <-raster(raster_file)# Create coordinates for all grid cellscoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")coords$PUID <-1:nrow(coords)#---------------------- FISHING DATA PROCESSING ----------------------## Convert aggregated_data to an sf objectfishing_sf <- aggregated_data %>%filter(!is.na(lon_05deg), !is.na(lat_05deg), !is.na(predicted_fishing_hours)) %>%st_as_sf(coords =c("lon_05deg", "lat_05deg"), crs =4326)# Make sure CRS matchesst_crs(fishing_sf) <-st_crs(meow_ecos)# Process fishing data for ECOREGIONSfishing_with_ecoregion <-st_join(fishing_sf, meow_ecos %>% dplyr::select(ECOREGION, REALM))ecoregion_fishing_stats <- fishing_with_ecoregion %>%st_drop_geometry() %>%group_by(ECOREGION, REALM) %>%summarize(mean_fishing_hours =mean(predicted_fishing_hours, na.rm =TRUE),median_fishing_hours =median(predicted_fishing_hours, na.rm =TRUE),fishing_q05 =quantile(predicted_fishing_hours, 0.05, na.rm =TRUE),fishing_q95 =quantile(predicted_fishing_hours, 0.95, na.rm =TRUE),fishing_cells =n(),.groups ='drop' ) %>%filter(!is.na(ECOREGION)) %>%arrange(desc(mean_fishing_hours))#---------------------- PRIORITY DATA PROCESSING ----------------------## Combine all priority data with coordinatescombined_data <- CAPTAIN2_EDGE2_data %>%rename(Priority_EDGE2 = Priority) %>%left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, Priority) %>%rename(Priority_FUSE = Priority), by ="PUID") %>%left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, Priority) %>%rename(Priority_IUCN = Priority), by ="PUID") %>%left_join(coords, by ="PUID")# Convert to sf object for each indexiucn_sf <- combined_data %>%filter(!is.na(Longitude), !is.na(Latitude), !is.na(Priority_IUCN)) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =4326)fuse_sf <- combined_data %>%filter(!is.na(Longitude), !is.na(Latitude), !is.na(Priority_FUSE)) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =4326)edge2_sf <- combined_data %>%filter(!is.na(Longitude), !is.na(Latitude), !is.na(Priority_EDGE2)) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =4326)# Make sure CRS matchesst_crs(iucn_sf) <-st_crs(meow_ecos)st_crs(fuse_sf) <-st_crs(meow_ecos)st_crs(edge2_sf) <-st_crs(meow_ecos)# Process priority data for ECOREGIONS - IUCNiucn_with_ecoregion <-st_join(iucn_sf, meow_ecos %>% dplyr::select(ECOREGION, REALM))ecoregion_iucn_stats <- iucn_with_ecoregion %>%st_drop_geometry() %>%group_by(ECOREGION, REALM) %>%summarize(mean_priority =mean(Priority_IUCN, na.rm =TRUE),median_priority =median(Priority_IUCN, na.rm =TRUE),priority_q05 =quantile(Priority_IUCN, 0.05, na.rm =TRUE),priority_q95 =quantile(Priority_IUCN, 0.95, na.rm =TRUE),priority_cells =n(),.groups ='drop' ) %>%filter(!is.na(ECOREGION)) %>%arrange(desc(mean_priority))# Process priority data for ECOREGIONS - FUSEfuse_with_ecoregion <-st_join(fuse_sf, meow_ecos %>% dplyr::select(ECOREGION, REALM))ecoregion_fuse_stats <- fuse_with_ecoregion %>%st_drop_geometry() %>%group_by(ECOREGION, REALM) %>%summarize(mean_priority =mean(Priority_FUSE, na.rm =TRUE),median_priority =median(Priority_FUSE, na.rm =TRUE),priority_q05 =quantile(Priority_FUSE, 0.05, na.rm =TRUE),priority_q95 =quantile(Priority_FUSE, 0.95, na.rm =TRUE),priority_cells =n(),.groups ='drop' ) %>%filter(!is.na(ECOREGION)) %>%arrange(desc(mean_priority))# Process priority data for ECOREGIONS - EDGE2edge2_with_ecoregion <-st_join(edge2_sf, meow_ecos %>% dplyr::select(ECOREGION, REALM))ecoregion_edge2_stats <- edge2_with_ecoregion %>%st_drop_geometry() %>%group_by(ECOREGION, REALM) %>%summarize(mean_priority =mean(Priority_EDGE2, na.rm =TRUE),median_priority =median(Priority_EDGE2, na.rm =TRUE),priority_q05 =quantile(Priority_EDGE2, 0.05, na.rm =TRUE),priority_q95 =quantile(Priority_EDGE2, 0.95, na.rm =TRUE),priority_cells =n(),.groups ='drop' ) %>%filter(!is.na(ECOREGION)) %>%arrange(desc(mean_priority))#---------------------- MERGE DATASETS FOR SCATTERPLOTS ----------------------## Merge ECOREGION data for each indexecoregion_combined_iucn <-inner_join( ecoregion_fishing_stats, ecoregion_iucn_stats, by =c("ECOREGION", "REALM")) %>%mutate(cell_ratio = fishing_cells / priority_cells)ecoregion_combined_fuse <-inner_join( ecoregion_fishing_stats, ecoregion_fuse_stats, by =c("ECOREGION", "REALM")) %>%mutate(cell_ratio = fishing_cells / priority_cells)ecoregion_combined_edge2 <-inner_join( ecoregion_fishing_stats, ecoregion_edge2_stats, by =c("ECOREGION", "REALM")) %>%mutate(cell_ratio = fishing_cells / priority_cells)#---------------------- CREATE ECOREGION SCATTERPLOTS ----------------------## IUCN - Get top 5 fishing and top 5 priority ecoregionstop_fishing_iucn <- ecoregion_combined_iucn %>%arrange(desc(mean_fishing_hours)) %>%head(5) %>%pull(ECOREGION)top_priority_iucn <- ecoregion_combined_iucn %>%arrange(desc(mean_priority)) %>%head(5) %>%pull(ECOREGION)ecoregions_to_label_iucn <-unique(c(top_fishing_iucn, top_priority_iucn))label_data_iucn <- ecoregion_combined_iucn %>%filter(ECOREGION %in% ecoregions_to_label_iucn)# Create IUCN plotiucn_ecoregion_plot <-ggplot(ecoregion_combined_iucn,aes(x = mean_fishing_hours, y = mean_priority)) +geom_point(aes(size = fishing_cells, color = REALM), alpha =0.4) +geom_label_repel(data = label_data_iucn,aes(label = ECOREGION),size =3,max.overlaps =10,box.padding =0.5 ) +scale_x_log10() +scale_size_continuous(name ="Number of Cells", range =c(1, 8)) +scale_color_manual(values =c("red", "blue", "darkgreen", "purple", "orange", "brown", "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +labs(title ="Relationship Between Fishing Pressure and Conservation Priority",subtitle ="By Ecoregion (IUCN)",x ="Mean Fishing Hours (log scale)",y ="Mean Conservation Priority (IUCN)") +theme_minimal() +guides(color =guide_legend(title ="Realm", override.aes =list(size =4)))# FUSE - Get top 5 fishing and top 5 priority ecoregionstop_fishing_fuse <- ecoregion_combined_fuse %>%arrange(desc(mean_fishing_hours)) %>%head(5) %>%pull(ECOREGION)top_priority_fuse <- ecoregion_combined_fuse %>%arrange(desc(mean_priority)) %>%head(5) %>%pull(ECOREGION)ecoregions_to_label_fuse <-unique(c(top_fishing_fuse, top_priority_fuse))label_data_fuse <- ecoregion_combined_fuse %>%filter(ECOREGION %in% ecoregions_to_label_fuse)# Create FUSE plotfuse_ecoregion_plot <-ggplot(ecoregion_combined_fuse,aes(x = mean_fishing_hours, y = mean_priority)) +geom_point(aes(size = fishing_cells, color = REALM), alpha =0.4) +geom_label_repel(data = label_data_fuse,aes(label = ECOREGION),size =3,max.overlaps =10,box.padding =0.5 ) +scale_x_log10() +scale_size_continuous(name ="Number of Cells", range =c(1, 8)) +scale_color_manual(values =c("red", "blue", "darkgreen", "purple", "orange", "brown", "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +labs(title ="Relationship Between Fishing Pressure and Conservation Priority",subtitle ="By Ecoregion (FUSE)",x ="Mean Fishing Hours (log scale)",y ="Mean Conservation Priority (FUSE)") +theme_minimal() +guides(color =guide_legend(title ="Realm", override.aes =list(size =4)))# EDGE2 - Get top 5 fishing and top 5 priority ecoregionstop_fishing_edge2 <- ecoregion_combined_edge2 %>%arrange(desc(mean_fishing_hours)) %>%head(5) %>%pull(ECOREGION)top_priority_edge2 <- ecoregion_combined_edge2 %>%arrange(desc(mean_priority)) %>%head(5) %>%pull(ECOREGION)ecoregions_to_label_edge2 <-unique(c(top_fishing_edge2, top_priority_edge2))label_data_edge2 <- ecoregion_combined_edge2 %>%filter(ECOREGION %in% ecoregions_to_label_edge2)# Create EDGE2 plotedge2_ecoregion_plot <-ggplot(ecoregion_combined_edge2,aes(x = mean_fishing_hours, y = mean_priority)) +geom_point(aes(size = fishing_cells, color = REALM), alpha =0.4) +geom_label_repel(data = label_data_edge2,aes(label = ECOREGION),size =3,max.overlaps =10,box.padding =0.5 ) +scale_x_log10() +scale_size_continuous(name ="Number of Cells", range =c(1, 8)) +scale_color_manual(values =c("red", "blue", "darkgreen", "purple", "orange", "brown", "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +labs(title ="Relationship Between Fishing Pressure and Conservation Priority",subtitle ="By Ecoregion (EDGE2)",x ="Mean Fishing Hours (log scale)",y ="Mean Conservation Priority (EDGE2)") +theme_minimal() +guides(color =guide_legend(title ="Realm", override.aes =list(size =4)))# Create the plots without titles and legends (we'll use a shared legend)# IUCN Plotiucn_ecoregion_plot <-ggplot(ecoregion_combined_iucn,aes(x = mean_fishing_hours, y = mean_priority)) +geom_point(aes(size = fishing_cells, color = REALM), alpha =0.4) +geom_label_repel(data = label_data_iucn,aes(label = ECOREGION),size =3,max.overlaps =10,box.padding =0.5 ) +scale_x_log10() +scale_size_continuous(name ="Number of Cells", range =c(1, 8)) +scale_color_manual(values =c("red", "blue", "darkgreen", "purple", "orange", "brown", "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +labs(x ="Mean Fishing Hours (log scale)",y ="Mean Conservation Priority") +theme_minimal() +theme(legend.position ="none", # Remove legend from individual plotspanel.grid.minor =element_blank(),plot.background =element_rect(fill ="white", color =NA),panel.background =element_rect(fill ="white", color =NA) ) +guides(color =guide_legend(title ="Realm", override.aes =list(size =4, alpha =1)),size =guide_legend(title ="Number of Cells") )# FUSE Plotfuse_ecoregion_plot <-ggplot(ecoregion_combined_fuse,aes(x = mean_fishing_hours, y = mean_priority)) +geom_point(aes(size = fishing_cells, color = REALM), alpha =0.4) +geom_label_repel(data = label_data_fuse,aes(label = ECOREGION),size =3,max.overlaps =10,box.padding =0.5 ) +scale_x_log10() +scale_size_continuous(name ="Number of Cells", range =c(1, 8)) +scale_color_manual(values =c("red", "blue", "darkgreen", "purple", "orange", "brown", "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +labs(x ="Mean Fishing Hours (log scale)",y ="Mean Conservation Priority") +theme_minimal() +theme(legend.position ="none", # Remove legend from individual plotspanel.grid.minor =element_blank(),plot.background =element_rect(fill ="white", color =NA),panel.background =element_rect(fill ="white", color =NA) ) +guides(color =guide_legend(title ="Realm", override.aes =list(size =4, alpha =1)),size =guide_legend(title ="Number of Cells") )# EDGE2 Plotedge2_ecoregion_plot <-ggplot(ecoregion_combined_edge2,aes(x = mean_fishing_hours, y = mean_priority)) +geom_point(aes(size = fishing_cells, color = REALM), alpha =0.4) +geom_label_repel(data = label_data_edge2,aes(label = ECOREGION),size =3,max.overlaps =10,box.padding =0.5 ) +scale_x_log10() +scale_size_continuous(name ="Number of Cells", range =c(1, 8)) +scale_color_manual(values =c("red", "blue", "darkgreen", "purple", "orange", "brown", "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +labs(x ="Mean Fishing Hours (log scale)",y ="Mean Conservation Priority") +theme_minimal() +theme(legend.position ="none", # Remove legend from individual plotspanel.grid.minor =element_blank(),plot.background =element_rect(fill ="white", color =NA),panel.background =element_rect(fill ="white", color =NA) ) +guides(color =guide_legend(title ="Realm", override.aes =list(size =4, alpha =1)),size =guide_legend(title ="Number of Cells") )# Create one plot with legend for extracting the common legendlegend_plot <-ggplot(ecoregion_combined_iucn,aes(x = mean_fishing_hours, y = mean_priority)) +geom_point(aes(size = fishing_cells, color = REALM), alpha =0.4) +scale_size_continuous(name ="Number of Cells", range =c(1, 8)) +scale_color_manual(values =c("red", "blue", "darkgreen", "purple", "orange", "brown", "black", "pink", "darkgray", "navy", "darkred", "forestgreen")) +guides(color =guide_legend(title ="Realm", override.aes =list(size =4, alpha =1)),size =guide_legend(title ="Number of Cells") ) +theme_minimal() +theme(plot.background =element_rect(fill ="white", color =NA),panel.background =element_rect(fill ="white", color =NA) )# Create a blank plot for the legend positionblank_plot <-ggplot() +theme_void()# Extract the legend as a separate groblegend_grob <- cowplot::get_legend(legend_plot)# Use ggpubr::ggarrange to combine plots in a 2x2 grid with legend in position 4library(ggpubr)combined_ecoregion_plots <-ggarrange( iucn_ecoregion_plot, fuse_ecoregion_plot, edge2_ecoregion_plot, as_ggplot(legend_grob),labels =c("(A)", "(B)", "(C)", ""), # No label for legend positionncol =2, nrow =2)# Ensure white background for the entire combined plotcombined_ecoregion_plots <- combined_ecoregion_plots +theme(plot.background =element_rect(fill ="white", color =NA))# Print the combined plotprint(combined_ecoregion_plots)
Code
# Save the plot with adjusted dimensions for 2x2 layout and explicit white backgroundggsave(here::here("outputs", "ecoregion_fishing_priority_relationship.png"), combined_ecoregion_plots, width =12, height =10, dpi =300, bg ="white")# Print summary statisticscat("Summary of ecoregion-level analysis:\n")
# Optional: Create individual plots if needed# print(iucn_ecoregion_plot)# print(fuse_ecoregion_plot) # print(edge2_ecoregion_plot)# The main combined plot is printed above
Maps of SR, FUn and EDGE2
Code
# Sepcies richness ----# Load required librarieslibrary(here)library(terra)library(ggplot2)library(viridis)library(sf)library(rnaturalearth)library(smoothr)# Set the path to your raster filesraster_path <-here("Data", "tif files continental")# Get list of all TIF files in the directorytif_files <-list.files(raster_path, pattern ="\\.tif$", full.names =TRUE)# Check if files were foundif(length(tif_files) ==0) {stop("No TIF files found in the specified directory")}print(paste("Found", length(tif_files), "TIF files"))
[1] "Found 1000 TIF files"
Code
# Read all rasters into a SpatRaster stackspecies_rasters <-rast(tif_files)# Calculate species richness by summing across all layers# This assumes presence = 1, absence = 0 or NAspecies_richness <-sum(species_rasters, na.rm =TRUE)# Name the layernames(species_richness) <-"Species_Richness"# Create a basic plotplot(species_richness, main ="Species Richness Map",col =viridis(100))
Code
# Alternative: Create a ggplot map with McBryde-Thomas 2 projection# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Project the species richness raster to McBryde-Thomas 2species_richness_projected <-project(species_richness, mcbryde_thomas_2)# Convert to dataframe for ggplotrichness_df <-as.data.frame(species_richness_projected, xy =TRUE)# Remove NA values for cleaner plottingrichness_df <- richness_df[!is.na(richness_df$Species_Richness), ]# Project the world mapworld_projected <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() )# Convert richness data to sf points for proper clippingrichness_sf <-st_as_sf(richness_df, coords =c("x", "y"), crs = mcbryde_thomas_2)# Clip the richness data to the globe boundaryrichness_clipped <-st_intersection(richness_sf, globe_border)# Extract coordinates back for plottingrichness_clipped_coords <-st_coordinates(richness_clipped)richness_final <-data.frame(x = richness_clipped_coords[,1],y = richness_clipped_coords[,2],Species_Richness = richness_clipped$Species_Richness)# Create ggplot mapspecies_richness_plot <-ggplot() +geom_tile(data = richness_final, aes(x = x, y = y, fill = Species_Richness)) +geom_sf(data = world_projected, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_fill_viridis_c(name ="Species\nRichness",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +coord_sf(crs = mcbryde_thomas_2, expand =FALSE) +labs(title ="Global Species Richness Map",x =NULL, y =NULL) + my_theme# Display the plotprint(species_richness_plot)
Code
# Richness of threatened species, top 25% FUSE and top 25% EDGE2 ----# Load required librarieslibrary(here)library(terra)library(ggplot2)library(viridis)library(sf)library(rnaturalearth)library(smoothr)library(dplyr)# Load the conservation metrics datasetconservation_data <-read.csv(here("Data", "My dataframes", "continental_shark_conservation_metrics_10_harmonised_IUCN_csv.csv"))# Filter for threatened species (IUCN values of 0.5, 0.75, or 1)threatened_species <- conservation_data %>%filter(IUCN %in%c(0.5, 0.75, 1)) %>%pull(Species.name)# Get top 25% FUSE scoring species by rank (not by threshold)n_species <-nrow(conservation_data)top_25_percent_count <-ceiling(n_species *0.25)top_fuse_species <- conservation_data %>%arrange(desc(FUSE)) %>%slice_head(n = top_25_percent_count) %>%pull(Species.name)# Get top 25% EDGE2 scoring species by ranktop_edge2_species <- conservation_data %>%arrange(desc(EDGE2)) %>%slice_head(n = top_25_percent_count) %>%pull(Species.name)# Get the actual threshold values for reportingfuse_threshold <- conservation_data %>%arrange(desc(FUSE)) %>%slice(top_25_percent_count) %>%pull(FUSE)edge2_threshold <- conservation_data %>%arrange(desc(EDGE2)) %>%slice(top_25_percent_count) %>%pull(EDGE2)print(paste("Total species in dataset:", n_species))
print(paste("Found", length(top_fuse_species), "top 25% FUSE species (threshold:", round(fuse_threshold, 4), ")"))
[1] "Found 250 top 25% FUSE species (threshold: 0.1 )"
Code
print(paste("Found", length(top_edge2_species), "top 25% EDGE2 species (threshold:", round(edge2_threshold, 4), ")"))
[1] "Found 250 top 25% EDGE2 species (threshold: 0.1 )"
Code
# Set the path to your raster filesraster_path <-here("Data", "tif files continental")# Get list of all TIF files in the directoryall_tif_files <-list.files(raster_path, pattern ="\\.tif$", full.names =TRUE)# Extract species names from file names (assuming filename format includes species name)# You may need to adjust this pattern based on your actual file naming conventionfile_species_names <-basename(all_tif_files)file_species_names <-gsub("\\.tif$", "", file_species_names)# Function to find matching raster files for a species listfind_matching_files <-function(species_list, all_files) { matching_files <-c()for(species in species_list) {# Create pattern to match species name in filename (handles spaces and underscores) species_pattern <-gsub(" ", "[_ ]", species) matches <- all_files[grepl(species_pattern, basename(all_files), ignore.case =TRUE)] matching_files <-c(matching_files, matches) }return(unique(matching_files))}# Find matching files for each species groupthreatened_files <-find_matching_files(threatened_species, all_tif_files)fuse_files <-find_matching_files(top_fuse_species, all_tif_files)edge2_files <-find_matching_files(top_edge2_species, all_tif_files)print(paste("Found", length(threatened_files), "raster files for threatened species"))
[1] "Found 366 raster files for threatened species"
Code
print(paste("Found", length(fuse_files), "raster files for top FUSE species"))
[1] "Found 250 raster files for top FUSE species"
Code
print(paste("Found", length(edge2_files), "raster files for top EDGE2 species"))
[1] "Found 250 raster files for top EDGE2 species"
Code
# Check if files were foundif(length(threatened_files) ==0) {stop("No matching TIF files found for threatened species. Please check file naming convention.")}# Read rasters and calculate richness for each group# Threatened speciesthreatened_rasters <-rast(threatened_files)threatened_richness <-sum(threatened_rasters, na.rm =TRUE)names(threatened_richness) <-"Threatened_Species_Richness"# Top FUSE speciesfuse_rasters <-rast(fuse_files)fuse_richness <-sum(fuse_rasters, na.rm =TRUE)names(fuse_richness) <-"Top_FUSE_Richness"# Top EDGE2 speciesedge2_rasters <-rast(edge2_files)edge2_richness <-sum(edge2_rasters, na.rm =TRUE)names(edge2_richness) <-"Top_EDGE2_Richness"# Create basic plots for all three metricsplot(threatened_richness, main ="Threatened Species Richness Map",col =viridis(100))
Code
plot(fuse_richness, main ="Top 25% FUSE Species Richness Map",col =viridis(100))
Code
plot(edge2_richness, main ="Top 25% EDGE2 Species Richness Map",col =viridis(100))
Code
# Function to create projected richness mapcreate_richness_map <-function(richness_raster, title, subtitle, legend_name) {# Project the richness raster to McBryde-Thomas 2 richness_projected <-project(richness_raster, mcbryde_thomas_2)# Convert to dataframe for ggplot richness_df <-as.data.frame(richness_projected, xy =TRUE)colnames(richness_df)[3] <-"richness_value"# Remove NA values for cleaner plotting richness_df <- richness_df[!is.na(richness_df$richness_value), ]# Convert richness data to sf points for proper clipping richness_sf <-st_as_sf(richness_df, coords =c("x", "y"), crs = mcbryde_thomas_2)# Clip the richness data to the globe boundary richness_clipped <-st_intersection(richness_sf, globe_border)# Extract coordinates back for plotting richness_clipped_coords <-st_coordinates(richness_clipped) richness_final <-data.frame(x = richness_clipped_coords[,1],y = richness_clipped_coords[,2],richness_value = richness_clipped$richness_value )# Create ggplot map richness_plot <-ggplot() +geom_tile(data = richness_final, aes(x = x, y = y, fill = richness_value)) +geom_sf(data = world_projected, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_fill_viridis_c(name = legend_name,guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +coord_sf(crs = mcbryde_thomas_2, expand =FALSE) +labs(title = title,subtitle = subtitle,x =NULL, y =NULL) + my_themereturn(richness_plot)}# Create maps for all three metricsthreatened_richness_plot <-create_richness_map( threatened_richness, "Global Threatened Species Richness Map","Species with IUCN threat levels 0.5, 0.75, or 1.0","Threatened\nSpecies Richness")fuse_richness_plot <-create_richness_map( fuse_richness,"Global Top 25% FUSE Species Richness Map", paste("Species in top quartile of FUSE scores (≥", round(fuse_threshold, 4), ")"),"Top FUSE\nSpecies Richness")edge2_richness_plot <-create_richness_map( edge2_richness,"Global Top 25% EDGE2 Species Richness Map",paste("Species in top quartile of EDGE2 scores (≥", round(edge2_threshold, 4), ")"),"Top EDGE2\nSpecies Richness")# Display all plotsprint(threatened_richness_plot)
Code
print(fuse_richness_plot)
Code
print(edge2_richness_plot)
Differences between priority and normal richness maps
Code
create_comparison_maps <-function(captain2_data, richness_raster, priority_threshold =0.85, comparison_name ="Comparison") {# Count grid cells with CAPTAIN2 priority > threshold high_priority_cells <-sum(captain2_data$Priority > priority_threshold, na.rm =TRUE)print(paste("CAPTAIN2", comparison_name, "cells with priority >", priority_threshold, ":", high_priority_cells))# Get the same number of top richness cells richness_values <-values(richness_raster, na.rm =TRUE) richness_threshold <-quantile(richness_values, 1- (high_priority_cells /length(richness_values)), na.rm =TRUE)print(paste("Richness threshold for top", high_priority_cells, "cells:", round(richness_threshold, 2)))# Create binary rasters for comparison# CAPTAIN2 binary (1 = high priority, 0 = other) captain2_binary <- captain2_data captain2_binary$Binary_Priority <-ifelse(captain2_data$Priority > priority_threshold, 1, 0)# Richness binary (1 = top cells, 0 = other) richness_binary <- richness_rastervalues(richness_binary) <-ifelse(values(richness_raster) >= richness_threshold, 1, 0)# Convert CAPTAIN2 data to sf if it's not alreadyif(!inherits(captain2_binary, "sf")) {# Assuming captain2_data has Longitude and Latitude columns captain2_binary_sf <-st_as_sf(captain2_binary, coords =c("Longitude", "Latitude"), crs =4326) } else { captain2_binary_sf <- captain2_binary }# Transform to match the richness raster CRS captain2_binary_sf <-st_transform(captain2_binary_sf, crs =crs(richness_binary))# Convert to raster using terra::rasterize with explicit method captain2_raster <-rasterize(captain2_binary_sf, richness_binary, field ="Binary_Priority", fun ="max", # Use max to avoid averagingbackground =0) # Set background to 0# Ensure we only have 0s and 1s captain2_raster[captain2_raster >0] <-1# Define the McBryde-Thomas 2 projection mcbryde_thomas_2 <-"+proj=mbt_s"# Calculate difference and debug each step difference_raster <- captain2_raster - richness_binaryprint(paste("Difference raster total cells:", ncell(difference_raster)))print(paste("Difference raster non-NA cells:", sum(!is.na(values(difference_raster)))))# Create a mask for cells where at least one method has priority (value = 1) priority_mask <- (captain2_raster ==1) | (richness_binary ==1)print(paste("Priority mask TRUE cells:", sum(values(priority_mask), na.rm =TRUE)))# Apply the mask to keep only priority cells difference_raster[!priority_mask] <-NAprint(paste("After masking, non-NA cells:", sum(!is.na(values(difference_raster)))))# Project difference raster difference_projected <-project(difference_raster, mcbryde_thomas_2)print(paste("After projection, non-NA cells:", sum(!is.na(values(difference_projected)))))# Convert to dataframe and clip to globe diff_df <-as.data.frame(difference_projected, xy =TRUE)colnames(diff_df)[3] <-"difference" diff_df <- diff_df[!is.na(diff_df$difference), ]print(paste("After df conversion, rows:", nrow(diff_df)))# Convert to sf and clip diff_sf <-st_as_sf(diff_df, coords =c("x", "y"), crs = mcbryde_thomas_2)print(paste("SF object rows:", nrow(diff_sf)))# Create world map and globe border for this function world <-ne_countries(scale ="medium", returnclass ="sf") world_projected <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding box and border globe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90)) globe_border <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base theme my_theme <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank() ) diff_clipped <-st_intersection(diff_sf, globe_border)print(paste("After globe clipping, rows:", nrow(diff_clipped)))# Extract coordinates diff_coords <-st_coordinates(diff_clipped) diff_final <-data.frame(x = diff_coords[,1],y = diff_coords[,2],difference = diff_clipped$difference )print(paste("Final data rows:", nrow(diff_final)))# Create difference map diff_plot <-ggplot() +geom_tile(data = diff_final, aes(x = x, y = y, fill =factor(difference))) +geom_sf(data = world_projected, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_fill_manual(values =c("-1"="red", "0"="black", "1"="blue"),labels =c("-1"="Only Richness", "0"="Agreement", "1"="Only CAPTAIN2"),name ="Priority\nDifference",guide =guide_legend(title.position ="top", title.hjust =0.5) ) +coord_sf(crs = mcbryde_thomas_2, expand =FALSE) +labs(title =paste("Priority Comparison:", comparison_name),subtitle =paste("CAPTAIN2 vs Simple Richness (", high_priority_cells, "top cells each)"),x =NULL, y =NULL) + my_theme +theme(legend.position ="bottom")# Calculate summary statistics total_cells <-nrow(diff_final) agreement_cells <-sum(diff_final$difference ==0, na.rm =TRUE) captain2_only <-sum(diff_final$difference ==1, na.rm =TRUE) richness_only <-sum(diff_final$difference ==-1, na.rm =TRUE) agreement_percent <- (agreement_cells / total_cells) *100 captain2_only_percent <- (captain2_only / total_cells) *100 richness_only_percent <- (richness_only / total_cells) *100# Print summary statistics with debugging infocat("\n========================================\n")cat("SUMMARY STATISTICS -", comparison_name, "\n")cat("========================================\n")cat("CAPTAIN2 priority cells:", high_priority_cells, "\n")cat("Richness priority cells:", high_priority_cells, "(same number)\n")cat("CAPTAIN2 raster 1s:", sum(values(captain2_raster) ==1, na.rm =TRUE), "\n")cat("Richness raster 1s:", sum(values(richness_binary) ==1, na.rm =TRUE), "\n")cat("Total priority cells compared:", total_cells, "\n")cat("----------------------------------------\n")cat("Agreement (both prioritize same cell):", agreement_cells, "(", round(agreement_percent, 2), "%)\n")cat("Only CAPTAIN2 prioritizes:", captain2_only, "(", round(captain2_only_percent, 2), "%)\n")cat("Only Richness prioritizes:", richness_only, "(", round(richness_only_percent, 2), "%)\n")cat("----------------------------------------\n")cat("Total accounted for:", round(agreement_percent + captain2_only_percent + richness_only_percent, 2), "%\n")cat("========================================\n\n")return(list(plot = diff_plot, captain2_cells = high_priority_cells,richness_threshold = richness_threshold,total_cells = total_cells,agreement_cells = agreement_cells,agreement_percent = agreement_percent,captain2_only_percent = captain2_only_percent,richness_only_percent = richness_only_percent))}# Threatened species comparisonthreatened_comparison <-create_comparison_maps( CAPTAIN2_IUCN_data_nonzero, threatened_richness, priority_threshold =0.85,comparison_name ="Threatened Species")